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ABSTRACT 

O ; Structure formation in the intergalactic medium (IGM) produces large-scale, 

; collisionless shock waves, where electrons can be accelerated to highly relativis- 

j> ; tic energies. Such electrons can Compton scatter cosmic microwave background 

photons up to 7-ray energies. We study the radiation emitted in this process 
using a hydrodynamic cosmological simulation of a ACDM universe. The re- 
sulting radiation, extending beyond TeV energies, has roughly constant energy 
flux per decade in photon energy, in agreement with the predictions of Loeb & 
00 1 Waxman (2000). Assuming that a fraction £ e = 0.05 of the shock thermal energy 

is transferred to the population of accelerated relativistic electrons, as inferred 
from collisionless non-relativistic shocks in the interstellar medium, we find that 

O 



in 



the energy flux of this radiation, e 2 (dJ/de) ~ 50 — 160 eV cm~ 2 s^ 1 sr _1 , consti- 
tutes ~ 10% of the extragalactic 7-ray background flux. The associated 7-ray 
point-sources are too faint to account for the ~ 60 unidentified EGRET 7-ray 
sources, but GLAST should detect and resolve several 7-ray sources associated 



Of 



2 ' with large-scale IGM structures for £ e ~ 0.03, and many more sources for larger 

£ e . The intergalactic origin of the shock-induced radiation can be verified through 
a cross-correlation with, e.g., the galaxy distribution that traces the same struc- 
ture. Its shock-origin may be tested by cross-correlating its properties with radio 
synchrotron radiation, emitted as the same accelerated electrons gyrate in post- 
shock magnetic fields. We predict that GLAST and Cherenkov telescopes such 
as MAGIC, VERITAS and HESS should resolve 7-rays from nearby (redshifts 
z < 0.01) rich galaxy clusters, perhaps in the form of a ~ 5 — 10 Mpc diameter 
ring-like emission tracing the cluster accretion shock, with luminous peaks at its 
intersections with galaxy filaments detectable even at z ~ 0.025. 
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1. Introduction 



In the past three decades, clear evidence has emerged indicating the existence of an 
extragalactic 7-ray background (EGRB). The origin of this radiation, however, has remained 
highly speculative. The first unambiguous detection of isotropic extragalactic 7-ray emission 
was obtained by the SAS-2 satellite at 1977. Subsequent experiments, especially the EGRET 
instrument aboard the Compton Gamma Ray Observatory (CGRO), have confirmed the 
existence of this radiation. These measurements indicate a generally isotropic (fluctuation 
amplitude < 20% on > 20° scales) spectrum in the energy range 30 MeV-120 GeV, well fitted 
by a single power law, with photon number density per energy interval dJ/ale ~ e -( 21 o±o.03) 
(Sreekumar et al. 1998). Known extragalactic sources, such as blazars, account for less 
than 7% of the EGRB, and unidentified blazars probably contribute no more than 25% of 
the radiation (Mukherjee & Chiang 2000; but for a different model see Stecker & Salamon 
2001). Thus, most of the EGRB flux originates from a source yet to be identified. 

Recently, Loeb & Waxman (2000) have proposed a model which attributes most of the 
EGRB flux to emission from large-scale structure-forming regions in the universe. Strong 
collisionless, non-relativistic shock waves characteristic of these regions accelerate particles 
to highly relativistic energies by the Fermi-acceleration mechanism, first discussed as sources 
of ultra-high energy cosmic rays by Norman, Melrose & Achterberg (1995) and by Kang, 
Ryu & Jones (1996). Electrons, accelerated in such shocks to highly relativistic energies with 
Lorentz factors up to 7 max ~ 10 7 , scatter a small fraction of the cosmic microwave background 
(CMB) photons up to 7-ray energies, as inverse- Compton radiation. The estimated energy 
flux from this process per decade in photon energy, e 2 (dJ/de) « 1.5 keV cur 2 s _1 sr _1 , is in 
good agreement with the detected EGRB, for typical cosmological parameters, provided that 
the shocked matter accumulated a present-day mass-average temperature T ~ 10 7 K, and 
that a fraction £ e « 0.05 of the associated thermal energy was transferred to the relativistic 
electrons. In this model, nearby rich galaxy clusters are expected to be bright 7-ray sources, 
as they host strong, dense shocks (Loeb & Waxman 2000; Totani & Kitayama 2000). 

Kawasaki & Totani (2001) report a correlation between a sample of high-latitude steady 
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unidentified EGRET sources and pairs or groups of galaxy clusters, where strong accretion 
and merger shocks are expected. However, this sample is too small and the correlation too 
weak for it to be statistically significant. EnBlin et al. (2001) suggest that various features 
of the asymmetric radio plasma ejected from the giant radio galaxy NGC 315 can provide 
the first indication of a structure-formation shock in a filament of galaxies. According to the 
above model, such a shock should be accompanied by 7-ray emission, yet to be detected. 

The model presented above has several interesting implications, other than an identifi- 
cation of the origin of the EGRB: 7-ray emission could be the first tracer of the undetected 
warm-hot intergalactic medium (WHIM), gas at temperatures 10 5 K < T < 10 7 K which is 
believed to contain a large fraction of the baryons in the low redshift universe (e.g. Dave 
et al. 2001a). The angular fluctuations in the EGRB should, according to the model, re- 
flect the underlying forming structure, and thus serve as a diagnostic of intergalactic shocks. 
The relativistic electrons emit synchrotron radiation, as they gyrate in intergalactic mag- 
netic fields. The resulting radio radiation, although weaker than the CMB, could have a 
fluctuating component dominating the CMB fluctuations at low (v < 10 GHz) frequencies 
(Waxman & Loeb 2000). Inverse- Compton emission from shock-accelerated electrons in the 
hard X-ray (HXR) and extreme ultra-violet (EUV) bands was recently studied by Miniati et 
al. (2001). They find that such emission could account for a small fraction of the measured 
EUV excess, and for the entire HXR excess in a few reported galaxy clusters, provided that 
the latter contain relatively strong merger or accretion shocks. 

In order to test this model, calibrate its free parameters and extract useful predictions, 
one must follow the evolution of structure in the intergalactic medium (IGM) in detail. The 
complex, non-linear processes and many scales involved, make this a difficult problem, ren- 
dering numerical simulations a powerful and indispensable tool. In this paper, we apply 
cosmological simulations to the problem of non-thermal emission from shock-accelerated 
electrons, and demonstrate how various EGRB predictions can be efficiently extracted from 
them. For this purpose, we have developed the necessary tools to localize intergalactic shocks 
in the simulations, inject the corresponding shock-accelerated electrons into the gas, and cal- 
culate the resulting radiation. We consider a hydrodynamic simulation of a ACDM universe, 
currently considered the most successful theory of structure formation, which describes a flat 
universe containing dark matter and presently dominated by vacuum energy. The relativis- 
tic electrons are assumed to carry a fraction £ e of the shock thermal energy. We adopt the 
value £ e = 0.05, which we show is inferred from collisionless non-relativistic shocks in the 
interstellar medium. 

Overall, we find an emitted inverse- Compton spectrum well-described by a mildly bro- 
ken power law, dJ/de ~ e - 2 - 04 i n the energy range 10 keV — 1 GeV and dJ/de ~ e ~ 213 
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in the energy range 10 GeV — 5 TeV, in good agreement with the prediction of Loeb & 
Waxman (2000). The calculated energy flux (per decade in photon energy), e 2 (dJ/de) = 
50 — 160 eV cm" 2 s _1 sr -1 , constitutes ~ 10% of the EGRB flux and is smaller by a factor 
of ~ 10 than estimated above, mainly due to the lower present-day temperature reached by 
the simulation, To ~ 4 x 10 6 K, and its slower heating rate. We predict the existence of 7-ray 
halos, associated with large-scale structure. Although too faint to account for the EGRET 
unidentified 7-ray sources, several of these halos should be detectable by future experiments, 
such as the Gamma-Ray Large-Area Space Telescope (GLAST, planned 4 to be launched in 
2006), provided that £ e > 0.03. The number of detectable sources is sensitive to the fraction 
of shock-energy transferred to the relativistic electrons, with ~ 15 well-resolved GLAST 
sources predicted for £ e = 0.05. Detection of intergalactic shock-induced 7-ray sources may 
be used to calibrate the free parameter £ e . We examine images of such an object, as predicted 
to be produced by GLAST and by the Major Atmospheric Gamma-ray Imaging Cherenkov 
(MAGIC) telescope (expected 5 to become operational during 2002). These images reveal 
an elliptic ring-like feature, of diameter corresponding to 5 — 10 Mpc, surrounding a nearby 
[z ~ 0.012) rich galaxy cluster. This ring, tracing the cluster accretion shock, has luminous 
peaks along its circumference, probably indicating the locations of intersections with galaxy 
filaments, channeling gas into the cluster. 

We describe the Loeb- Waxman model in §2. We first consider the gross characteristics of 
intergalactic shocks based on order of magnitude estimates. In this context, we parameterize 
the distribution of accelerated electrons, and discuss the inverse- Compton and synchrotron 
radiation they emit. The simulation we chose to study and its basic predictions are described 
in §3. We discuss the underlying theory and cosmological model, present some character- 
istics of the simulated universe and examine its structure. Section 4 describes how EGRB 
predictions can be extracted from cosmological simulations in general and then focuses on 
the method applied to the simulation discussed in §3. In §5, we present our results regard- 
ing the radiation predicted from the simulation. We describe the resulting inverse-Compton 
spectrum, analyze the corresponding 7-ray sky maps and study the features and distribution 
of the simulated 7-ray sources. In §6 we discuss the implications of our results for future 
experiments and models. Appendix A contains some formulae used to calculate and inte- 
grate the radiation emitted by the accelerated electrons. Appendix B presents an algorithm 
devised to identify point sources in our simulated maps of the 7-ray sky. 



4 See http://glast.gsfc.nasa.gov 

5 See http://hegral.mppmu.mpg.de/MAGICWeb 



- 5 - 



2. Theory 



This section describes the model for non-thermal emission from structure formation 
shocks proposed by Loeb and Waxman (2000), using simple order-of-magnitude estimates. 
First, we estimate the parameters of the large-scale shock waves involved in structure forma- 
tion. Next, we discuss the high-energy electrons accelerated by such shocks. We evaluate the 
power index and cutoff of their distribution, and evaluate its normalization by estimating 
the efficiency of electron acceleration in intergalactic shocks. We conclude by calculating the 
radiation emitted by these electrons, mainly through inverse-Compton scattering off CMB 
photons. The estimates presented in this section needed for our numerical simulation are 
discussed in more detail in §4. 



2.1. Intergalactic Shock Waves 



Large-scale structure in the universe is believed to have evolved by the gravitational 
collapse of initially over-dense regions. As such an initial seed accreted increasing amounts 
of mass, large converging flows of matter were accelerated towards it, inevitably brought to 
a violent halt. This process resulted in large-scale shock waves, forming as massive gas flows 
collided with opposite flows or with the newly formed object. In the following, we concentrate 
our attention on the 'recent' stages of structure formation, taking place at moderate to low 
redshifts (z < 3), and deduce some order of magnitude estimates of the resulting shock wave 
properties, following Cen & Ostriker (1999). 

Linear theory of structure formation predicts that the length and mass scales of an 
object, forming (becoming non-linear) at these low redshifts, are roughly given by: 

\nl(z) = A and M NL (z) = M f{zf , (1) 

1 + z 

where Ao and M are the length and mass scales of structures forming at the present time 
and f(z) is a slowly varying function of the redshift z, satisfying f(z — 0) = 1, such that 
smaller scales become non-linear first. These parameters are sensitive to the cosmological 
model. We use a 'concordance' ACDM model of Ostriker & Steinhardt (1995) - a flat universe 
with normalized vacuum energy density ft\ = 0.7, matter energy density Qm = 0.3, baryon 
energy density Qb = 0.04, Hubble parameter h = 0.67, and initial perturbation spectrum of 
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slope n = 1 and normalization <t 8 = 0.9 - to find: 

A ~ lOh^ 1 Mpc and M ~ 1.5 x lO 14 /^ 1 (jf^j M Q , (2) 

where h 67 = fa/0. 67. The order of magnitude estimates obtained from linear theory in 
this subsection should be regarded with caution, as non-linear effects become important in 
the forming objects, and the phases of the perturbations can no longer be ignored. The 
distribution of the forming objects according to linear theory may be better estimated using 
the approach of Press & Schechter (1974). 

Geometric considerations dictate that the gas in an object of size A, collapsing over a 
time period t colu will achieve velocities of order v ~ \/t coU . Since a structure that became 
non-linear at redshift z will collapse over a time period of order t con ~ H{z)~ 1 ) where H(z) 
is the Hubble parameter, the formation of such an object is expected to involve shocks with 
velocity: 

v sh (z) ~ k X NL (z) H{z) = v MMfl > ( 3 ) 

where k is a dimensionless number of order unity, vq = kXoHo is the typical shock velocity 
at z = and h(z) = H(z)/H(z = 0). For the ACDM model described above we find: 

v ~ 700 km s" 1 and h(z) = ^tt A + f2 M (l + zf . (4) 

Hence, the shock waves resulting from converging flows during structure formation in the 
IGM have a non-relativistic velocity similar to those of shocks surrounding supernova rem- 
nants (where accelerated particles are observed) and they are steadily growing in both scale 
and velocity. 

The baryonic component may be treated as an ideal gas with an adiabatic index Y = 
5/3, and is sufficiently rarefied such that the shocks concerned are collisionless. Using the 
Rankine-Hugoniot shock adiabat (Landau & Lifshitz 1959), we can relate the shocked gas 
temperature to the shock velocity, 



T(z) 



f(z)h(z) 



l 2 



(5) 



1 + z 

where T = [m p / (Tk B )](av ) 2 is the temperature of gas, shocked at z — 0, ks is the Boltz- 
mann constant, and a = c Sj d/v sh is the ratio between the post-shock (downstream) sound 
velocity and the shock velocity (relative to the unshocked gas), approaching ~ 0.56 for a 
strong shock. For the ACDM model above we find: 

T ~ 10 7 K . (6) 

As we show in §3.2, this estimate of T is indeed valid only to order of magnitude and appears 
to exceed the true value by a factor ~ 3. 
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2.2. Shock- Accelerated Electron Population 



2.2.1. Acceleration Mechanism 

Collisionless, non-relativistic shock waves are generally known to accelerate a power- 
law energy distribution of high-energy particles. This phenomenon has been observed in 
astrophysical shock waves on various scales, such as in shocks forming when the supersonic 
solar wind collides with planetary magnetospheres, in shocks surrounding supernovae (SNe) 
remnants in the interstellar medium, and probably also in shocks in many of the most active 
extragalactic sources: quasars and radio galaxies (Drury 1983; Blandford & Eichler 1987). 
These power law distributions extend up to energies ~ 100 keV in the Earth's bow shock, 
to ~ 100 TeV in SNe remnants (Tanimori et al. 1998) where shock velocities are similar to 
those of intergalactic shocks, and over 10 6 TeV in Galactic cosmic rays. 

Magnetic fields of amplitude B ~ 0.1 — 1 /iG are measured in halos of galaxy clusters 
(Kronberg 1994) and in the Coma super-cluster (Fusco-Femiano et al. 1999), suggesting that 
the IGM in which the shocks propagate is significantly magnetized. The measured magnetic 
field energy constitutes a fraction £ B — 0.01 of the typical post-shock thermal energy in 
cluster environments (Waxman & Loeb 2000). 

Electrons are expected to be accelerated in intergalactic shocks by the Fermi acceleration 
mechanism. The velocity difference between the two sides of a shock front acts as a first- 
order acceleration mechanism idEjdt ~ E), provided that slow- moving elastic scattering 
centers are embedded in the gas. Alfven waves or magnetic inhomogeneities can function 
as such centers, scattering electrons through small angles and causing them to repeatedly 
diffuse through a shock front. Diffusion of the energetic electrons downstream, away from 
the shock front, provides a Poissonian escape route. These effects lead to a characteristic 
power-law distribution of the high-energy electrons steady-state number density: 

<hlp (E) ~ AE~ S , (7) 



alE, 



e 



where A is a normalization constant. The vast size and perpetual duration of these shocks 
enable them to dissipate significant fractions of their energies in this process. 

The resulting distribution of high energy electrons can be calculated using the test- 
particle approximation (Drury 1983; Blandford & Eichler 1987), where the test particles, 
accelerated by the shock, are assumed to have no effect on its structure. This approximation, 
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valid for the acceleration of a small population of high energy electrons, predicts that the 
power-law index of the resulting population depends only on the kinematic structure of the 
shock front, and is given by: 

where r is the compression factor of the shock. For a shock in an ideal gas, r is given by: 

r ~ f+T + (r + 1)M 2 ' ^ ' 

where M = v sh /c StU is the Mach number, and v sh and c s>u denote the shock velocity and the 
pre-shock (upstream) sound velocity, respectively. The shock waves of interest are generally 
strong, i.e. involve high (M ^> 1) Mach numbers. Hence, we find r ~ (r + l)/(r — 1) = 4, 
and thus s xs 2. 

In order to completely define the power-law distribution discussed, we must specify its 
normalization and the energy range over which it extends. The Fermi acceleration mechanism 
per-se does not impose stringent limits on the electron energies achieved, as long as the 
total energy carried by the accelerated electrons is small compared to the shock energy and 
the shock structure remains intact. The e-folding time for electron acceleration is given 
(Loeb & Waxman 2000) by: 

r«~T~2x lOS^-T^j)" 1 yr , (10) 

where r L is the electron Larmor radius, 77 is the electron Lorentz factor - 7 e - in units of 
10 7 , B_ 7 is the magnetic field in units of 0.1 /xG, and T d 7 is the downstream temperature 
of the gas, in 10 7 K. This time-scale is much shorter than the shock lifetime. The limit on 
electron energy is therefore not determined by the acceleration process, but rather by cooling 
processes: electrons will stop accelerating once their cooling rate overcomes their acceleration 
rate. The dominant cooling process for the relevant parameters is inverse-Compton scattering 
of the electrons off CMB photons, with a characteristic cooling time: 

r IC = - 3meC ~ 2 x lOV^l + z)-' Yr , (11) 
4cr T 7 e w CMB 

where u CMB is the CMB energy density. Synchrotron cooling of the electrons is less efficient, 
with a characteristic cooling time r syn = (u cmb /u b )t ic ~ 10 3 (£>_ 7 )~ 2 (1 + z) A r lC) where u B is 
the magnetic field energy density. Thus, inverse-Compton cooling is faster, as long as the 
downstream magnetic field satisfies B < 3(1 + z) 2 /iG. Equating r acc and r IC leads to the 
maximal Lorentz factor: 

7max = 3.3 x 10 7 (B_ 7 T dJ ) 1/2 (1 + z)- 2 . (12) 
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Since the electron energy loss is dominated by scattering of CMB photons, the flux of 
high energy photons produced by shock accelerated electrons depends only weakly on the 
magnetic field strength. The magnetic field determines, nevertheless, the highest energy 
of accelerated electrons, and hence the energy to which the inverse- Compt on spectrum ex- 
tends. For B ~ 0.1 /iG, as observed in cluster halos and which corresponds to a fraction 
£ B ~ 0.01 of the post-shock thermal energy, the highest energy electrons up-scatter CMB 
photons to energy e ~ 1 TeV. It should be pointed out here, that the observed sub-/xG 
fields may be representative of only the downstream (post shock) magnetic field strength, 
and that the upstream field may be significantly lower. In this case, the electron acceleration 
time will be longer, and the highest energy of accelerated electrons and of inverse- Compt on 
photons may be significantly lower than estimated from equation (12). The strength of the 
upstream magnetic field depends on the processes leading to IGM magnetization, which are 
not yet known. Magnetic fields may be produced by turbulence induced by the large scale 
structure shocks themselves (Kulsrud et al. 1997), by "contamination" by the first gener- 
ation of stars (Rees 1987) or radio sources (Daly & Loeb 1990; Furlanetto & Loeb 2001), 
or by galactic winds (Kronberg, Lesch, & Hopp 1999). Under all these scenarios, the up- 
stream magnetic field strength is not expected to be much smaller than the downstream 
field strength. Moreover, near equipartition fields are expected to be produced in collision- 
less shocks through the growth of electromagnetic instabilities, as indicated, for example, by 
observed supernovae remnants (Helfand & Becker 1987; Cargill & Papadopoulos 1988) and 
7-ray bursts (Gruzinov & Waxman 1999; Medvedev & Loeb 1999), and there is evidence 
from supernovae remnant observations for an enhanced level of magnetic waves ahead of the 
shock front (Achterberg, Blandford, & Reynolds 1994) We therefore adopt for our estimates 
of the maximum electron energy magnetic field strengths corresponding to fixed ~ 0.01, 
which reproduce the observed large scale field strengths. 

We can now parameterize the normalization of the relativistic electron distribution, by 
assuming that a fraction £ e of the thermal energy density induced by the shock, Au th , is 
transferred to these electrons: 

3 

u e = £ e Au th = -iek B (n d T d - n u T u ) , (13) 

where u e = m e c 2 J^ ma -* r y(dn e /d'y) is the total energy density of relativistic electrons and 
u and d subscripts indicate upstream (pre-shock) and downstream (post-shock) values, cor- 
respondingly. The gas number density, n, is related to the baryon number density, n bar , 
by n(m) = n h!lI m P) where (m) ps 0.59m p is the average particle mass in an ionized gas 
with a hydrogen mass fraction x — 0.76. For strong shocks s ~ 2 and n u T u <C ndTd, thus 
u e ~ Aln7 max , and we find: 

3£ e n d k B T d 
2 In 7 max 
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2.2.2. Acceleration Efficiency 

The efficiency of particle acceleration in structure formation shocks, parameterized 
above by £ e , is an important parameter bearing linearly on our results, and thus deserves a 
special discussion. Although no existing model credibly calculates the acceleration efficiency, 
we can evaluate £ e in intergalactic shocks by estimating the acceleration efficiency of other 
astrophysical shock waves. 

The best analogy to intergalactic shocks may be found in supernova remnant (hereafter 
SNR) shocks, drawing on the similarity between the velocities of intergalactic shocks and 
SNR shocks, both of the order of ~ 10 3 km s" 1 . The physics of shock waves is essentially 
determined by three parameters: the shock velocity v sh , the pre-shock density n u and the 
pre-shock magnetic field B u . Although the plasma density is very different in intergalactic 
shocks and in SNR shocks, the density may be extracted from the problem by measuring 
time in units of v~l, where v p ^ is the plasma frequency of the ions. The pre-shock mag- 
netic field, parameterized by the cyclotron frequency of the ions, z/ c i , may not be scaled 
out of the problem as well. However, by comparing v ci to the growth rate of electromag- 
netic instabilities in the shocked plasma, z/ om . = z/ Pi jf sh /c, we find that for strong shocks 
(u c i/z/e.m.) 2 = jz ,3^ 87r i <C I. We thus assume there is a well behaved limit when this ratio 
approaches zero, suggesting that the pre-shock magnetic field has little effect on the charac- 
teristics of strong shocks. With this assumption, we expect to find much similarity between 
strong shocks in different environments, as long as their shock velocities are comparable. 

The energy of the relativistic electrons accelerated in SNR shocks was calculated by 
several authors, using the multi-frequency emission from SNR's, preferably remnants with 
TeV detection such as SN1006 (G327.6 - 1.4; Tanimori et al. 1998). Dyer et al. (2001) have 
modeled the multi-frequency emission from SN1006, finding a shock efficiency of £ e = 5.3%. 
The total energy they find in relativistic electrons constitutes a fraction r] e ~ 1.4% out of 
the total supernova explosion energy (W ~ 5 x 10 50 erg, as often quoted in the literature). 
Other authors (Mastichiadis & de Jager 1996; Aharonian & Atoyan 1999) have estimated 
r] e ~ 1% — 2% in SN1006, which corresponds to £ e values in the range of 3.8% — 7.6%. Ellison, 
Slane and Gaensler (2001) have modeled the emission from supernova remnant G347.3 — 0.5, 
estimating the energy of relativistic electrons to constitute at least 1.2% - and more likely 
2.5% - of the shock kinetic energy flux, corresponding to 2.5% < £ e ~ 5%. 

An independent, less reliable method for estimating the electron acceleration efficiency 
relies on the observation, that relativistic ions accelerated in SNR shocks are dynamically 
significant (Reynolds & Ellison 1992; Baring et al. 1999), and thus carry a large fraction 
(tens of percents, e.g. ~ 50% according to Ellison et al. 2001) of the post shock thermal 
energy. The efficiency of electron acceleration may then be found by estimating the ratio r] 
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between the relativistic electron energy and the energy of relativistic ions. One may easily 
estimate the value of 77 in interstellar medium cosmic rays: using the electron-to-proton ratio 
in a given energy, e.g. 1% at GeV energies (Barwick et al. 1998), correcting for electron 
cooling and integrating over energies, we find rj ~ 3% — 4%, i.e. £ e ~ 1% — 2%. However, 
since there is no clear relation between 77 in the interstellar medium and its value immediately 
behind the ~ 10 3 km s _1 shocks found in SNR's and in the intergalactic medium, this 
estimate of £ e is controversial. 

To summarize, a shock efficiency of £ e ~ 0.05 is found from studies of supernovae 
remnants, characterized by shock velocities similar to those of intergalactic shocks and thus 
expected to exhibit similar behavior. Since there is some uncertainty in this parameter, 
values of £ e in the range 0.02 — 0.10 are considered plausible. 

2.3. Resulting Radiation 



In this subsection we calculate the radiation emitted by the relativistic electron pop- 
ulation described above. As mentioned earlier, most of the energy is produced when the 
electrons inverse-Compton scatter CMB photons, resulting in a spectrum that extends up 
to TeV energies. A secondary radiation process is synchrotron emission from the electrons, 
gyrating in the magnetic field, resulting in a similarly sloped spectrum extending up to the 
infrared. Other radiative processes experienced by the electrons, such as Bremsstrahlung 
and Cherenkov radiation, are less important. 

We begin by calculating the diffuse radiation, resulting from the various intergalactic 
shocks at recent epochs (moderate to low redshifts). We then focus our attention on a single 
forming structure and the angular distribution of relevant sources. 

2.3.1. Diffuse Radiation 



Inverse-Compton scattering describes the interaction of a free electron with radiation, 
when the electron kinetic energy exceeds the photon energy. In the ultra-relativistic limit, 
7 e ^> 1, the average relative increase in photon energy is given by Ae 7 /e 7 ~ (4/3)7^. Hence, 
an electron with energy E e = r y e m e c 2 ^> m e c 2 will scatter a CMB photon to an average 
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energy: 

<e 7 ) ~ ^ 2 e e CMB = CE 2 e , (15) 

where we have defined ( = (4/3) (m e c 2 ) 2 e CMB , and e CMB — 3.83kBT CMB is the average energy 
of a CMB photon. The extent of the spectrum is determined by the electron range of 
energies: The maximal electron Lorentz factor, 7 max ~ 3.3 x 10 7 , is given by equation (12), 
and equation (11) imposes a minimal energy threshold, as electrons with Lorentz factors 
smaller than 7 min ~ 200 hardly radiate during a Hubble time ~ 10 10 yr. Plugging this range 
of Lorentz factors into equation (15) shows that the emitted photon spectrum extends from 
~ 50 eV up to ~ TeV energies. 

The emitted spectrum may be found by assuming that each relativistic electron scatters 
a photon only once. This approximation, discussed in §4.3, gives qualitatively correct results 
and is exact for strong shocks. With this simplification, we find that the emitted energy 
density per unit photon energy is given by: 

d u i lr \ _ 1 dn e ( /—f-:\ _ A- ts j2-\ -sji 



"(0 = 7^ VVC )= i?' i - l e-« i . (16) 



de 7 w 2(dE e V v ' V 2 

For a strong shock, where the electron number density satisfies dn e (E)/dE e ~ AE~ 2 , this 
reduces to: 

e^-(e)= A/2. (17) 

For weak shocks, the approximation in equation (16) introduces a small numerical error 
(~ 2), because the temporal evolution of the electron distribution affects the normalization 
of the resulting spectrum. 

Determining the present-day photon spectrum requires integration over different shocks 
at various times. This, in turn, is determined by the redshift-dependence of the IGM pa- 
rameters, and requires the use of non-linear structure formation theory. We note that the 
contribution of a shock that developed at redshift z to the observed value of e[du 1 (e)/de 1 ], 
is roughly proportional to (1 + z)~ 4 . Thus, the main contribution to the spectrum observed 
today, especially at high photon energies, is expected to originate from recent shocks at 
z < 0.5. Assuming that a fraction / sh of the baryons in the universe were recently shocked 
to a temperature T ~ 10 7 K, when 7 max already approached 10 7 , we find from equation (14) 
and equation (17): 

e ^( e ) ^ 6 .4 x 10- 7 T 7 ( ( eV cm' 3 . (18) 
rfe 7 W 7 V0.05 J V 0.04 J v ; 

The inverse-Compton flux per decade in photon energy thus becomes: 

e 2 ^( e ) ~ 1.5 T 7 ( ( ^f) keV cm" 2 s" 1 sr" 1 . (19) 



rfe 7 w V 0.05 J V 0.04 
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This result is in excellent agreement with observational data, when supplemented by 
the contribution from unresolved point sources (Loeb & Waxman 2000). One can employ 
the Press-Schechter mass function in order to carry out the time-integration over the ra- 
diation resulting from intergalactic shocks at various epochs. This approach gives, for 
slightly different cosmological parameters (Qa = 0.65, Qb = 0.05, h = 0.7), a similar result 
(Loeb & Waxman 2000): 

d T IC ' PS f P \ 

e 2_)_( e ) ~ 1.5 \ —\ keV cm" 2 s" 1 sr" 1 . (20) 

The synchrotron radiation emitted by the accelerated electrons may now be estimated. 
A relativistic electron with a Lorentz factor 7 e , gyrating in a magnetic field of amplitude 
B, emits photons with characteristic energy e 7 ~ [(heB) / '(m e c)]7 2 . Synchrotron radiation is 
thus expected in photon energies ranging from 10 -12 eV and up to 1 eV. The emitted flux 
can be crudely estimated as: 

2 d J^ yn 2 d J^ c ub 



e 2 — Me) ~ e z — y -(e) 



- "^m)^)"™-"'- 1 "- 1 - (21) 

whereas a more careful approach is to assume that the energy density of the post-shock mag- 
netic field constitutes a fraction of the shock-induced thermal energy. This assumption, 
combined with the Press-Schechter mass function, leads to a similar (somewhat higher) flux: 

where a typical value £b = 0.01 produces magnetic fields consistent with observations of the 
outer halos of X-ray clusters (Waxman & Loeb 2000). 



2.3.2. Radiation From Forming Objects, Angular Statistics 



We start with the general case of an astrophysical object formed recently, and calculate 
its inverse- Compton luminosity. We assume that the object accretes mass at a rate M, 
heating it up to a temperature T by strong shocks. The baryon number accretion rate is thus 
given by M^Ib/^m)^)' 1 . We assume, as previously, that a fraction £ e of the post-shock 
thermal energy is converted into relativistic electrons. The maximal electron Lorentz factor 
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7max is determined from equation (12), and depends on the assumed magnetic field. Electrons 
with cooling time longer than the virialization time of the object, i vir , will not have sufficient 
time to radiate a significant fraction of their energy. This introduces a low-energy cutoff in 
the emitted spectrum, at photon energy e min = (4/3)7^ in e CMB , where 7 min ~2x 10 12 yr/t vir . 
This effect reduces the total luminosity by a factor In (7 max /7 min )/ hi7 max . Combining these 
considerations, we find that the total and specific luminosities of the object are given by: 

L 10 = M^(m)- 1 i e -k B T ln (7 m ax/7 min ) ? (23) 

2 m7 max 

and: 

^ = 2^£k3 fOT » 1 (W) 2 feV * * S (^ ™ . (24) 

Structure formation at the current epoch produces dense, large-scale sheets and fila- 
ments, entangled in space. The densest regions appear at the intersections of these features, 
around the locations of galaxy clusters. As an example, we consider such a galaxy cluster. 
The typical gas mass involved in the formation of a cluster is M gas ~ 1O 14 M , accreting 
over a virialization time t viT ~ 10 9 yr and shocked to a temperature T gas > 10 7 K. A mag- 
netic field of amplitude B ~ 0.1 //G, as suggested by observations, gives 7 min ~ 2000 and 
7 max ~ 3.3 x 10 7 . Whence, 



cluster 



and: 



22 - 1045 (4) 8w) («) fe) ~ erg »- • « 25) 

e 7 ; lustor ~ 0.05L^ uster for energies 5 keV < e 7 < 1 TeV . (26) 



de-y 



A comprehensive discussion of particle acceleration in galaxy clusters and the resulting 
inverse- Compton radiation up to hard X-rays, may be found in Sarazin (1999). 

The angular statistics of the radiation predicted by the model can be investigated by 
simplifying assumptions regarding the mass distribution. Waxman and Loeb (2000) cal- 
culated the angular correlation of the spectrum using the Press-Schechter mass function. 
On sub-degree angles they find fluctuations > 40% in the inverse- Compton spectrum and 
fluctuations > 100% in the synchrotron spectrum. 



3. Simulation 
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The previous section presented some order of magnitude estimates, leading from linear or 
non-linear structure formation theories to an EGRB spectrum, that roughly matches present- 
day observations. However, further progress along this route is severely compromised by the 
complexity and non-linearity of structure formation. This situation naturally calls for the 
use of cosmological simulations that incorporate hydrodynamical effects. For this project, we 
analyze a simulation of a ACDM cosmology. In the following, we describe the cosmological 
model and background theory of this simulation, and discuss some characteristics of the 
simulated universe and the predicted underlying structure. 



The cosmological model adopted in the simulation is the 'concordance' ACDM model 
of Ostriker & Steinhardt (1995), presented in §2.1. An initial Harrison-Zel'dovich spectrum 
of fluctuations was assumed, and linear structure formation theory (Eisenstein & Hu 1998) 
employed to predict the resulting spectrum at a certain redshift z atait , providing initial condi- 
tions for the simulation. The various parameters of the cosmological model are summarized 
in Table 1 (see e.g. Springel, White & Hernquist 2001b). 

The simulation was performed using the TreeSPH code GADGET (Springel, Yoshida, 
& White 2001a). This code, designed for both serial and parallel computers, is intended for 
collisionless and gas-dynamical cosmological simulations. 

Dark matter is modeled by the code as a self-gravitating collisionless fluid, represented by 
a large number N nM of particles. Newton's gravitational potential is modified by introducing 
a spline softening at small separations, equivalent to treating each particle as a normalized 
spline mass distribution, instead of as a point mass (see, e.g. Hernquist & Katz 1989). This 
modification is necessary in order to suppress large angle scattering in two-body collisions. 
Poisson's equation is solved using a hierarchical tree algorithm (Barnes & Hut 1986). 

The baryonic component of the universe is modeled by the simulation as an ideal gas, 
with an adiabatic index T = 5/3, and is represented by a large number of particles, similar 
to the dark matter. The gas is governed by the continuity equation: 



3.1. Model and Theory 




(27) 



the Euler (momentum) equation: 



p— = -VP - pV$ + V • E , 



(28) 
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and the thermal energy equation (the first law of thermodynamics): 

du P„ 1— „ A(u,p) , , 

= V • v - -VijdiVj - -^L , 29 

dt p p p 

where P, p, u and v are the pressure, mass density, thermal energy per unit mass and 
velocity of a gas element; $ is the gravitational potential, found from Poisson's equation: 

V 2 $(r,t) = 47rGp(r,t) ; (30) 

and we have used Lagrangian derivatives: 

T = # + v-V. (31) 
dt dt v 1 

An artificial viscosity (Steinmetz 1996) term E was added to the above equations, in order 
to capture shocks. A cooling function A(u, p) may be included in the energy equation, to 
incorporate radiative cooling processes. 

These hydro dynamical equations are solved using the SPH - Smoothed Particle Hydro- 
dynamics - technique (Lucy 1977; Gingold & Monaghan 1977; Monaghan 1992), well suited 
for non-axisymmetric three-dimensional astrophysical problems. The version of the code 
used for the simulation employs adaptive SPH smoothing lengths, such that each particle 
has a constant number N s of neighbors lying within its kernel. The hydrodynamical equa- 
tions and their solutions in the formalism of this version of SPH appear in Springel et al. 
(2001a). We note that the simulation analyzed in this paper employed a formulation of SPH 
which does not rigorously conserve entropy in adiabatic parts of the flow (for a discussion, 
see Hernquist 1993, Springel & Hernquist 2002a). 

The code simulates the evolution of dark matter and gas (hereafter: SPH) particles 
within a finite simulation box of comoving size £ 3 om , using periodic boundary conditions. 
Due to technical reasons, to be discussed in the next section, we have chosen to employ a 
simulation that did not include radiative cooling processes. Radiative cooling, dominated 
by Bremsstrahlung and line emission at the epoch of interest, is sub-dominant relative to 
inverse- Compton emission from relativistic electrons in all but the densest regions within 
the cores of clusters where the baryon number density exceeds ~ 10~ 3 cm -3 , and in cold 
(T < 10 4 K), dilute regions. We discuss the significance of radiative cooling towards the end 
of this section. The simulation neglected feedback from star formation and supernovae. The 
code parameters used in the simulation are presented in Table 2. 



3.2. Resulting Structure 
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The simulation chosen for the project, as most cosmological simulations, predicts a 
highly structured IGM in the present-day universe, as illustrated in Figure 1. Several hot 
and dense objects are spread throughout space, separated by large empty voids. The largest 
of these structures are sheet-like objects (pancakes) and filaments, entangled in space, some 
as long as ~ 100 Mpc. At the intersection of the filaments appear denser, spheroidal objects, 
with typical size ~ 1 Mpc, where galaxy clusters are formed. Simulations predict that the 
pancake-resembling objects are the first large structures to appear. At the intersections of 
these pancakes, filaments emerge, which in turn intersect and drain mass into spheroidal 
objects. At the large scales discussed, dark matter and baryons essentially trace the same 
structure. 

The objects described above are significantly denser than their surroundings. Thus, the 
last stages of their evolution are dominated by non-linear processes. These processes are 
mainly gravitational, though other physical processes, such as non-linear gas dynamics and 
radiative cooling (when included), influence these stages of structure evolution as well. 

It is instructive to examine the thermodynamic evolution of the universe, by applying 
some averaging schemes to the temperature and density at different epochs, as shown in 
Figure 2. The mass averaged temperature is notably higher than the volume averaged 
temperature, because hot regions are typically much denser, leaving most of the volume 
with a dilute, cooler gas (e.g. Dave et al. 1999). The present-day averaged temperature 
reached by the simulation, T(z = 0) = 3.9 x 10 6 K (see Figure 2), is significantly lower than 
the order-of-magnitude estimate presented earlier. (Note that the approximation of eq. [5] is 
in rough agreement with the mass averaged temperature, although deviating towards higher 
temperatures at recent times.) Furthermore, the simulation result for the present-day mass 
averaged temperature is lower than that predicted by Cen & Ostriker (1999) by a factor 
~ 2.5, but in good agreement with the results of e.g. Refregier et al. (2000), Croft et al. 
(2001), and Refregier & Teyssier (2000). R. Cen (private communication) informs us that 
revised estimates for the temperature evolution of the Cen & Ostriker (1999) simulation, 
correcting for an error in the integration of the thermal energy equation, now actually lie 
below the other simulation results. 

The volume averaged density is proportional to (1 + z) 3 , due to the expansion of the 
universe, whereas the approximately constant mass averaged density indicates that much 
of the mass has virialized. The temperature averaged density is higher than the volume 
averaged density, due to the above mentioned correlation between hot and dense regions. 
This (temperature averaged) density decreases in time, growing closer to the volume averaged 
density, as the dilute gas surrounding virialized structures heats up. 

We examine the effect of radiative cooling by super-imposing cooling-time contours on 
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the distribution of mass in the temperature-density phase-space (Figure 3, left). The figure 
suggests that radiative cooling has important effects on the present-day gas, residing in the 
most dense, hot regions and in the most cold, dilute regions. The cooling time in these 
volumes due to Bremsstrahlung and line emission, t cool , is comparable to or smaller than the 
dynamical time, t H ~ if -1 . Such regions, where t cool < t H , presently contain approximately 
25% of the gas mass, and a progressively larger mass fraction at higher redshifts (Figure 
3, right). The dense and hot regions, where radiative cooling is significant, are associated 
with the cores of clusters, well beyond the accretion shock fronts of interest and should thus 
have little effect on our results, although the weaker merger and flow shocks may be affected. 
Radiative effects in the cold, dilute regions cools the collapsing gas and may result in slightly 
stronger accretion shocks. However, as these regions are cold to begin with, this effect should 
be unimportant. 



4. Method 



Next we discuss how cosmological simulations can be exploited to yield various predic- 
tions regarding the EGRB, according to the Loeb-Waxman model. Parts of the discussion 
are specific to the cosmological simulation chosen for the project, introduced in the previous 
section, and to the EGRB extraction schemes used. Predictions of shock-induced radia- 
tion are obtained from the simulation in three basic steps: First, the simulation is used to 
keep track of the forming structures and locate the shocks involved in their virialization. 
Second, relativistic electron populations are injected into the shocked gas, according to the 
predictions of the Fermi acceleration mechanism. And third, the radiation emitted by the 
relativistic electrons is calculated and integrated to give measurable EGRB features. The 
organization of this section corresponds to these three steps. 



4.1. Extracting Shocks 



Since shock waves are discontinuities in the thermodynamic quantities, locating shocks 
involves identifying large gradients in the properties of the gas. Gradients of quantities such 
as pressure, temperature or velocity, can be calculated directly from a simulation. In our 
case, this requires SPH kernel interpolation, as described in the previous section. Next, 
one can search for gradients larger than some imposed thresholds. These would imply the 
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presence of shocks, stronger than a corresponding shock strength. The shock front could 
then be traced out, yielding the shock topology and the jump conditions across the shock, 
which can in turn be checked to see if they satisfy the Rankine-Hugoniot adiabat. 

Such an approach was adopted by Miniati et al. (2000), using simulations based on a 
version of the total variation diminishing scheme (TVD), designed to capture strong shocks 
(Ryu et al. 1993). They identified very strong accretion shocks with Mach numbers M of up 
to 10 3 around the filaments and halos of their simulations. These halos exhibit, in addition, 
weaker shocks (M ~ 3 — 10), originating from the merger of structures or from accretion 
shocks into old structures, already embedded within newly formed halos. Combined, these 
shocks lead to a complicated pattern of inter-winding shocks. Miniati et al. find that the 
jump conditions, imposed by a shock, are altered by adiabatic gravitational compression, 
especially when the shock is close to the core of a cluster. Since the limited simulation 
resolution renders these two effects inseparable, shock parameters were extracted under the 
assumption that gravitational effects are isothermal. 



4-1.1. Shock- Extraction Scheme 



Since we are interested more in the shocked gas than in the topology of shocks, we 
find SPH simulations sufficient for our needs, although they do provide limited shock front 
resolution. This enables us to adopt a simpler approach for locating shocks, appropriate for 
Lagrangian simulations such as SPH. This method involves following a gas element through- 
out the simulation, searching for rapid increases in its entropy. In cosmological simulations, 
entropy changes of the gas should result only from shock waves or from cooling processes. 
This is the reason we have chosen, for simplicity, to work with a simulation in which radia- 
tive cooling is neglected, shown in the previous section to yield reliable shocks. In such a 
simulation, the entropy of the gas should only increase, and do so in discrete steps, deci- 
sively indicating the presence of shocks. Ideally, one could then identify all entropy 'jumps' 
as shock waves, determining their space-time coordinates. The thermodynamic jump con- 
ditions across a shock could then be found, by comparing the thermodynamic parameters 
of the gas element before (Tj, p { ) and after (Tj, pf) passing through the shock, directly, or 
using the entropy-change, AS = C v In [(Tf/T i )(pf/p i ) < - 1 ~ rS) ~\ , induced by the shock. Using 
the Rankine-Hugoniot adiabat, we find that the compression factor r of a shock - related to 
its Mach number by equation (9) - can be determined from AS according to: 

AS, = C v l AS = In 



r(r + i)-(r-i) r 
(r + i)-r(r-i) . 



(32) 
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In SPH simulations, one can choose the SPH particles themselves as the gas elements 
to be followed, since they are of constant mass and the simulation is Lagrangian. This tech- 
nically simplifies our shock-locating scheme, as kernel interpolation becomes unnecessary. 
It is instructive to study the distribution of mass - i.e. SPH particles - that accumulated 
various entropy-changes during recent epochs. Figure 4 (left) gives this distribution in terms 
of the mass fraction f(ASl) = m~J:dm(ASl) /dAS* of SPH particles that accumulated a 
normalized entropy change AS'* between z = 2 and z — 0, where m tot is the total mass 
of SPH particles. First, note that a small fraction of the gas experienced negative entropy 
changes, most probably as a result of lack of strict entropy conservation in the absence of 
shocks, owing to the particular formulation of SPH used for this simulation (see, e.g. Springel 
& Hernquist 2002a). The entropy-difference distribution among such particles resembles a 
normal distribution, consistent with accumulating noise. The shape of the remaining distri- 
bution has an interesting structure, suggesting a pronounced population of 'weak' shocks, 
with M ~ 3 — 10, and a large population of stronger shocks, with Mach numbers up to a 
few 10 3 , in agreement with Miniati et al. (2000). 

In order to determine the location and temporal duration of shock waves, one must 
examine the simulated gas in sufficiently short time intervals. This is also necessary in order 
to identify multiple shocks suffered by an SPH particle, as the effect of consecutive shocks 
is not additive, although the entropy change they induce is. We chose to examine snapshots 
of the simulated gas, separated by the simulation-box light-crossing time, as illustrated in 
Appendix A. For the mass resolution of the simulation, the entropy-increase experienced by 
most shocked particles is spread over several of these snapshots, indicating that shorter time 
intervals are not required. With the time intervals that were chosen, the entropy-difference 
distribution between any two snapshots (Figure 4, right) is fairly constant throughout the 
period of interest. Based on this distribution, we impose a minimal threshold AS™ [n , such 
that only particles experiencing a larger entropy-increase between two consecutive snapshots 
are considered shocked. A threshold value AS*™ 1 ' 1 = 0.4 was found to be reasonable, con- 
sidering numerical noise, although it entails the danger of eliminating some of the mildly 
weak shocks (M < 10), if an SPH particle requires many snapshots to cross such a shock. 
The fact that only ~ 10% of the gas suffers entropy-changes larger than AS™ in simplifies the 
computational work considerably. Demanding that the entropy increases at least by n s AS™ in 
during n s = 2 consecutive snapshots, and dismissing shocked particles that did not reach 
a temperature T min = 10 5 K by the present time, further simplifies the analysis. Our final 
results do not strongly depend upon the choice of these threshold values near those specified 
above. Note that with finite resolution, SPH particles may pass through a shock in a period 
longer than the snapshot temporal resolution. The simulation must thus be traced down to 
z < values in order to capture all recent shocks. 
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The temperature jump across a shock determines the thermal energy (per particle) it 
adds to the gas, AU th /N = (T — 1)^ 1 A; B (T (J — T u ), and is thus vital for the calculation of 
shock-induced radiation. In principle, this jump could be found directly, by comparing the 
temperature of the shocked SPH particle before and after a shock takes place. However, 
with finite simulation resolution, this approach leads to errors, due to significant adiabatic 
heating or cooling of the particle in parallel to the shock. In some extreme cases, we even 
found the shocked particle cooler than it was prior to the shock. (A similar problem when 
radiative cooling is active has been identified by e.g. Hutchings & Thomas [2000] and Martel 
& Shapiro [2001].) A partial remedy for this problem is to assume that the entropy induced 
by the shock - and thus the compression factor - was correctly determined so that the state of 
the shocked particle can be deduced from its state prior to the shock, or vice versa. Indeed, 
neither the pre-shock state of the particle nor its post-shock state are accurately known, due 
to the inaccurate shock timing and ongoing adiabatic processes, leading to some inevitable 
error. We chose to regard the temperature of the gas, before the shock was detected, as the 
true pre-shock temperature: T u = Tj. Hence, we take: 

T d = Tir {T ~ l) e^ s * (33) 

as the post-shock temperature, where r and AS** are related by equation (32). This choice 
yields more realistic results than the alternative, T d = Tf, as discussed in §4.1.2. 

4-1-2. Shock- Extraction Results 



The shock-extraction scheme, described above, identifies shocks in ~ 40% of the simu- 
lated gas mass during recent epochs [z < 2), where about 30% of the baryons were shocked 
at least once (Figure 5, left). Most identified shocks are concentrated around the halos in 
the simulation, although filaments and sheets are traced out as well. In some cases, but not 
always, shock fronts can be identified (Figure 6). 

The distribution of entropy-differences, accumulated by identified shocks (Figure 5, 
right), suggests that most of the strong shocks with Mach numbers in the range M > 10 2 
were located. It confirms, however, our suspicion of loosing most of the 'weak' shock popula- 
tion (M ~ 3 — 10). In such shocks, the entropy-change suffered by a single particle is spread 
over many snapshots, and is thus eliminated by the imposed entropy-change threshold. Fig- 
ure 5 (right) suggests that some of the very large accumulating entropy-increases, suffered 
by an SPH particle, result from multiple shocks. Separating the shocks in such cases results 
in a population of weaker shocks, artificially concentrated at the threshold shock strength 
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(M ~4). 

As a 'sanity check' of the proposed shock-location and energy-evaluation schemes, we 
compare the thermal energy gain of the simulation to the evaluated energy, induced by the 
identified shocks. We find that the energy, injected by our shocks, constitutes a fraction 
/ sh ~ 2/3 of the total thermal energy gained by the simulated gas (see Table 3 and Figure 
5, left). For comparison, in the estimates of §2.1 we took f sh = 1. Miniati et al. (2000) find, 
for a slightly different ACDM model (A = 0.55, h = 0.6), that the energy passing through 
shocks with Mach numbers M > 10 2 between redshift z = 1.5 and today, approximately 
equals the present-day thermal energy of their simulation. 

The above results were obtained by estimating the energy of an SPH particle, gained 
when passing through a shock, based on its temperature before passing through the shock 
(i.e. taking T u = T,j). Estimating this energy gain based on the temperature after passing 
through the shock (i.e. taking Td = Tf) results in shock energies lower, on average, by a 
factor ~ 4.5. This result is a combined effect of the low resolution of the simulation and 
adiabatic cooling, taking place in parallel to the shock. 

We find the shock-extraction results described above reassuring, as they indicate that 
the proposed shock location-and-evaluation recipe captures most of the strong shocks, and 
measures their energy reasonably well. Better performance of this algorithm may be achieved 
by using simulations with a higher mass resolution and with a more accurate handling of adi- 
abatic flows (e.g. Springel & Hernquist 2002a). It should be stressed that the elimination of 
weak shocks by imposing an entropy threshold for shock identification has only a small effect 
on the predicted diffuse 7-ray radiation, especially at high energies. First, the consistency 
check described above verifies that the energy injected by weak shocks is small compared 
to the energy injected by strong shocks: accounting for all weak shocks may increase the 
integrated flux at most by a factor of ~ 1.5. Second, the contribution of weak shocks to 
the 7-ray emission, especially at high photon energies, is small due to the steep distribution 
of the electrons accelerated by such shocks, reducing their 7-ray efficiency. At 10 GeV, for 
example, shocks with M — 10 have 66% efficiency (compared to strong shocks), dropping 
below 50% for M < 8. Weak shocks could, nonetheless, increase the 7-ray brightness of 
individual nearby structures such as galaxy clusters experiencing merger events, mainly in 
low photon energies, thereby slightly altering the appearance of their inner regions, and as 
a result may somewhat enhance the predicted source number counts. 
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4.2. Accelerated Electron Population 



Here we rely on the arguments of §2.2, deriving expressions more general and more 
precise than the order of magnitude estimates presented earlier. These results are used 
to calculate the relativistic electron distributions produced by shocks, and their temporal 
evolution. Appropriate electron populations can thus be 'injected' into the shocked SPH 
particles to simulate radiative processes. 



4-2.1. Production 



We postulate that a fraction £ e of the post-shock density of the thermal energy induced 
by a shock, Au th = [l/(r — l)]ndkB(T d — T u ), is converted into a power law distribution of rel- 
ativistic electrons, dn e (E)/dE e = AE~ S . The power law index is related to the compression 
factor by s = (r + 2)/(r — 1), where T = 5/3 was used here and is assumed henceforth. Impos- 
ing a distribution cutoff at the maximal energy attainable by the electrons, -E max = 7 ma x^eC 2 , 
fixes the normalization: 

f 1/ In 7 max for s = 2 ; 

= x { {s _ 2 )(m e c*y-y [1 - 7 -fr 2) ] for s > 2 . (34) 

We estimate the maximal Lorentz factor of the electrons by equating their acceleration 
e-folding time to the cooling time due to inverse- Compton scattering, yielding: 

3(r-l) eB d / r + i 



7 max = o TZX-^Z k B T u + r- — -k B T d . (35) 

8 {m)c 2 a T u CMB \ 1-1 J 

The average mass per particle is given by 

A-Tfl 

< m> = Sx+lU. " °' 59m ' ' <36> 
where x is the hydrogen mass fraction, Xe is the ionization ratio (per nucleon) and the second 
equality holds for a fully ionized gas, where Xe = (l + x)/2, with the hydrogen mass fraction 
of the simulation, x — 0.76, assumed henceforth. 

We parameterize the post-shock magnetic field, B d , by assuming that a fraction of 
the post-shock thermal energy is transferred into the energy of this field, giving: 

1/2 



B d = 4.3 x 10 



-8 



0.01 J \10n av (z) 



[1 + zf' 2 Gauss , (37) 
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where n av (z) is the average particle number density at redshift z. Hence, 



7max = l.lxl0 7 (rT^ + T^/4) 1 / 2 



T, 



d,7 



n 



0.01/ \lOn av (z 

which for strong shocks, where r — > 4 and T M <C T^, reduces to: 

1/4 



1/4 



-5/4 



7 (strons) = 2.2 x 10 7 T, 3/ 7 4 

' max a. { 



0.01 



n 



10n.Jz) 



[1 + z) 



-5/4 



(38) 



(39) 



4.2.2. Evolution 



After the relativistic electrons diffuse away from the shock, they cool, mainly by inverse- 
Compton scattering off CMB photons. The electron distribution - initially a power law - 
may thus evolve into a different form, as energetic electrons cool faster: 

P m = ^ = -CE* , (40) 

where we have defined C = 4u cmb ctt / "im^c? . This leads to a partial differential equation for 
the electron distribution: 

-n E (E, t) = CE 2 —n E (E, t) + 2CEn E (E, t) , (41) 
at oE 

where the notation n E = dn e /dE e is introduced. We ignore adiabatic cooling due to the 
Hubble expansion since the cooling timescales relevant for our calculated 7-ray spectrum 
are all much shorter than the Hubble time. For the initial (power-law) condition we find the 
solution: 

= f AE- S (1 - CEty*-V for E < min {£ max , E cool (t)} ; 
El ' 1 10 for£>min{£ max ,£ cool (t)} , 1 ' 

where we have defined E cool (t) = 1/Ct. This solution implies that the electron distribution 
steepens (if s > 2) towards E cool (t), beyond which all electrons have already cooled off. For 
strong shocks, where s — > 2, the electron distribution is stationary for E < E cool (t), as equal 
numbers of electrons enter and leave each energy bin at any given time. 

The inverse-Compton cooling time of a relativistic electron is a function of its Lorentz 
factor 7, approximately given by t coo i(7) ~ 2.3 x 10 12 7 _1 (1 + z)~ 4 yr. This implies that 
electrons with 7 < 7 min (;z), where 7 min is calculated in equation (A7), do not have sufficient 
time to cool between redshift z and today, and will not significantly contribute to the emitted 
radiation. 
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4.3. Emitted Spectrum 

^.3. 1. Inverse Compton Emission 



In the previous subsection we have demonstrated how the electron distribution, car- 
ried by each simulated SPH particle, can be calculated at any given time. Thus, it is 
straightforward, if somewhat lengthy, to calculate the inverse- Compton spectrum, emitted 
by each SPH particle. The emitted energy density per unit photon energy is obtained by 
integrating the single particle emission function, over both electron and photon distributions 
(Rybicki & Lightman 1979): 



— — 

de 1 



e , t) = l„ TC f l K£o)rf£o / 7-^(7,0/ (j^) d; , (43) 



where z^(eo) is the specific number density of the incident photon field, assumed constant (in 
the time period of interest) and isotropic, f\x) is the emission function, given by: 

f(x) = 2xlnx + x + 1 -2x 2 , (44) 

and we have assumed Thomson scattering in the electron rest frame: 7eo <C m e c 2 . The 
incident CMB photon field is approximately blackbody radiation, whence: 

o 2/L3 3 

= 1 iur \ i and ^cmb^ 2.73(1 + ^)K. (45) 

exp {e /k B T CMB ) - 1 



The assumptions made above are all satisfied for the ultra-relativistic electrons (7 ^> 
100) of interest: The CMB photon field is isotropic and varies slowly in time, thus it can 
be treated as constant for the fast emission of such electrons. The Thomson condition is 
satisfied for 77 <C 57(1 + z) -1 , rendering Klein-Nishina corrections unnecessary. 

A straightforward integration of equation (43) may be replaced by a simpler approach. 
The typical lifetime of a shock, comparable to the Hubble time t sh ~ t H = H^ 1 = 1.5 x 
10 w h(z)^ 1 yr, is much longer than the cooling time of ultra-relativistic electrons. Hence, a 
population of such electrons, constantly produced by a well resolved shock, will resemble a 
steady source of radiation. This enables us to approximate the emission from these electrons 
as instantaneous, i.e. we assume that all ultra-relativistic electrons with 7 > 7 min (z) loose all 
their energy to radiation, immediately after passing through the shock front. The error thus 
introduced is small, for electrons satisfying £ coo i(7) < f S h ^ 7 > 155h(z)/(l + z) 4 . Further 
justification for this approximation may be found in the temporal resolution of the snapshots 
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used. The time between consecutive snapshots examined is determined by the simulation- 
box light-crossing time, of order t ros ~ 6.5 x 10 8 (1 + Z)" 1 yr. Thus, it is obviously pointless to 
follow the temporal evolution of electrons, satisfying t cool (7) < t ICS ^7 > 3.5x 10 3 (1 + z)~ 3 . 

In the 'instantaneous emission' approximation described above, we are interested only 
in the time-integrated specific power, emitted by the accelerated electrons: 

du\ c . . dP}° 



de 1 



(e)= —^(e,t)dt. (46) 

Jo ae ~i 



Since the only time-dependence in the emitted power of equation (43) appears in the distri- 
bution of electrons, dn e ( r y,t)/d'y e , we may change the order of integration and first calculate 
the time-integrated electron distribution: 

'■(E) / n E (E,t)dt^—^—E-^, (47) 



dE e y ' j tv ' ' c( s -iy 

where the second equality is obtained by taking the upper integration limit to infinity or 
to max{t sh , t cool (E)}, which is well justified for ultra-relativistic electrons. Integrating the 
power emitted by the evolving electron distribution, as in equation (43), may thus be re- 
placed by calculating the instantaneous emission from the integrated electron distribution of 
equation (47). 

Since the integrated electron distribution is a simple power law, we may use the well 
known result for inverse-Compton emission of a power-law electron distribution, scattering 
a blackbody photon field (Rybicki & Lightman 1979). This leads to: 



du ic r 

^-(e) = Q(s)A\k B T CMB / (m e c 2 ) 2 



s/2-1 



e -s/2 



= Q{s)Ae^ s/2 ~ 1 \l + zy/ 2 - 1 e- s l\ (48) 
where Q(s) is a numerical coefficient, defined by: 

^-g y (, + ^" + 2) r ( 3+ 0^( 3 + l)- « 49 » 

((s) is the Riemann zeta function, and we defined e IC = (m e c 2 )' / [k B T CMB (z = 0)] = 1.1 PeV. 
This approximate result holds for emitted photon energies e in the range 47^ in eo < e <C 
47ma X e o- Writing this as: 

3-6 {^) 2 keV « e « 3.2 (^^) 2 TeV , (50) 
indicates that the approximation holds for the entire relevant 7-ray regime. 
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It is interesting to compare the above results to the approximation used in §2.3.1, where 
only the first photon emitted by each electron was treated, thus neglecting the evolution 
of the electron distribution. The approximate result of equation (16) agrees with the last 
equality in equation (48), up to a function of the power index s. The two would become 
identical, if one was to replace Q(s) in the latter by 1/2, and take e IC = 0.22 PeV. For 
a strong shock, we indeed find that Q(s = 2) = 1/2, and precisely recover equation (17) 
- the early estimate for a strong shock. This is not surprising, recalling that the electron 
distribution for s = 2 is mostly stationary, rendering evolutionary effects insignificant. The 
numerical mismatch between equation (16) and equation (48) varies as function of s, but 
remains of order 2 for relevant values of s (s < 8). 



4-3.2. Space-Time Integration 



Two different space-time integration schemes were used. The first, line of sight inte- 
gration, calculates the specific intensity /(e) at a given direction (9, fa) on the sky, where 9 
is its angle with the z-axis and <fi is the angle of its projection on the x-y plane with the 
rr-axis. Such a scheme may be used to study the statistics of the predicted radiation, by cal- 
culating moments of the specific intensity: the mean intensity J(e) = (1(e)), the correlation 
of intensities at various angular separations, etc. The second scheme provides direction bin 
integration, calculating the flux predicted within a small region a = (9±, 9 2 \ <f>i, fa) on the 
sky, defined by Q\ < 9 < 9 2 and fa < < fa. The flux, related to the specific intensity by 
F(e) = f a 1(e) cos9 dQ, is required for sky maps and for source identification and number 
counts. 

Both integration schemes involve simulating an integration volume (IV, for short), back- 
tracking the propagation of radiation towards an imaginary observer. The IV starts at an 
arbitrarily chosen observer location today, and is moved away from the observer and back in 
time, at the speed of light. The radiation from sources (i.e. shocked SPH particles) within 
the IV is calculated and summed up, to yield the total radiation detected by the observer 
today. Since the simulation-box light-crossing time is much shorter than the Hubble time, 
periodic boundary conditions must be used to simulate the motion of the IV. 

We treat the crossing of a shock front by an SPH particle as a single source of radiation 
for our integration scheme and call it, for short, a shock. Let such a shock, with specific 
luminosity L sh (e) and redshift z, lie at coordinates (r, 9o, fa) of a comoving coordinate system 
centered upon the observer. The contribution of this shock to the specific flux F sh (e), detected 
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by the observer today, satisfies: 



where cLl{z) = rR (l + z) is the luminosity distance between source and observer, R is the 
present-day value of the cosmo logical scale factor, and we have defined e' = (1 + z)e. With 
the spatial and temporal resolution used, many (>> 10) shocks fall within the IV in each time 
step. This enables us to approximate the shocks as point sources, such that the contribution 
of a shock with space-time coordinates (r, t) to the detected flux is determined by: 

A[F(e,a)] = \ F ^ if < r '*) G IV ^ ; (52) 
1 otherwise. 

The specific intensity contributed by the shock is calculated by approximating the latter 
as a sphere of uniform specific brightness _B sh (e). The proper radius of this sphere is deter- 
mined by the mass density p(z) of the SPH particle, according to d pi (z) = [3M SPH /47rp(2;)] 1 ^ 3 , 
where M SPH is the mass of an SPH particle. The brightness of the sphere is related to its 
detected flux by: 



F„Je)=irB. 



AMI ' (53) 

where <1a{z) = rR /(l + z) is the angular diameter distance between the source and the 
observer. Equating this result and equation (51) gives: 

^ = ^M^w • (54) 

and the contribution of the source to the specific intensity at direction (9, 0) is given by: 

A [/( e 0)] = / B ^( e ) if the ra y ^) intersects the sphere; 
i. otherwise. 



We approximate the specific luminosity of a shock by: 

L Bh (e) ~ 



■ (56) 



where At sb = tf — ti is the period of time during which the shock was detected. Since our 
shock locating-and-timing scheme is limited by the temporal resolution of the snapshots, At sh 
tends to overshoot the true duration of the shock. This error is partially compensated by 
the integration schemes, as they are limited by the same snapshot resolution: although the 
approximate luminosity will tend to undershoot its true value, the probability of including 
the shocked particle in the IV is enhanced proportionally. 
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5. Results 



Next we present various features of the radiation, inverse-Compton emitted by the 
shock-accelerated electrons, according to the cosmological simulation studied. We begin by 
presenting the predicted spectrum and discussing its implications. Next, we present some 
maps of the simulated 7-ray sky and discuss the sources they are composed of. We conclude 
with source number-counts extracted from such maps. 



5.1. Predicted Spectrum 



The inverse-Compton spectrum, predicted by the simulation, extends from the optical 
band into the far 7-ray regime, with photon energies in the range 1 eV < e < 10 4 TeV 
(Figure 7). However, the high energy tail of this spectrum is modified by photon-photon 
pair-production interactions, rendering the universe non-transparent for e > 1 TeV photons. 
Photons of e > 100 TeV thus interact with CMB photons or, at much higher energies, 
with radio background photons. Lower energy (1 TeV < e < 100 TeV) photons interact 
mainly with the intergalactic infrared (IR) background. The e + e~ pairs, produced in such 
interactions, may inverse-Compton scatter CMB or IR photons, thus initiating a cascade 
process (Nikishov 1962; Gould & Schreder 1967; Protheroe & Stanev 1993). As a result, the 
> 1 TeV spectrum is significantly attenuated, whereas at lower photon energies, piling-up 
of the cascading photons slightly enhances (flattens) the spectrum. 

The resulting spectrum provides roughly equal energy flux per decade in photon energy, 
in the photon energy range 10 keV < e < 1 TeV. This flux ranges from e 2 rfJ 7 (e)/c?e 7 ~ 
160 eV cm" 2 s" 1 sr" 1 at e = 10 keV, through 100 eV cm" 2 s" 1 sr" 1 at e = 1 GeV, and 
down to 50 eV cm" 2 s" 1 sr" 1 at e = 1 TeV. The number flux of photons in the energy range 
10 keV < e < 1 GeV is well approximated by: 

whereas in the energy range 10 GeV < e < 5 TeV the spectrum is slightly steeper and may 
be approximated, before including the cascade process, by: 

,iJ l( e ) ~ 8.8 x 10" 19 ( 6 ) ^ (eV cm 2 s sr) _1 . (58) 
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The slope of the resulting spectrum and the energy range over which it extends are in good 
agreement with the order of magnitude estimates, carried out in §2.3.1. The resulting flux, 
on the other hand, falls short of the early estimates of equation (19) and equation (20), by 
a factor ~ 10 (Figure 7). 

The discrepancy between the above flux, calculated from the simulation, and the flux 
predicted in §2.3.1, is due to a combination of several factors. First, the temperature reached 
by the simulation is a factor of ~ 2.6 lower than the order-of-magnitude temperature 10 7 K 
assumed previously. The thermal energy gain in the epoch examined (0 < z < 2) is thus 
smaller than the thermal energy used in §2.3.1 by a factor of ~ 3. Second, the energy 
induced by the identified shocks constitutes a fraction / sh ~ 2/3 of the thermal energy gain 
in the simulation, as opposed to the previous assumption / sh = 1. Third, the simulation 
identifies only approximately half of the shock-induced electron-energy in electrons, that 
had sufficient time to cool (7 > 7 min ), in contrast to the factor 1 — In (7 min )/ In (7 max ) ~ 2/3 
of §2.3.1, indicating the presence of weak shocks and very recent shocks. Finally, in §2.3.1 
we neglected the redshift of the emitted energy, caused by the expansion of the universe, 
assuming that most the detected radiation was emitted recently (z < 0.5). However, the 
energy accumulated by the detected shocks increases almost steadily for z < 2 (Figure 5, 
left), such that a significant fraction of the predicted flux originates at high redshifts: only 
~ 40% of this flux was emitted at z < 0.5. This effect reduces the predicted flux by a factor 
~ 5/3. Summing up the effects listed above, found from Table 3, yields the discrepancy 
factor: 3 • (3/2) • (4/3) • (5/3) = 10. 



5.2. Sky Maps, 7-ray Sources and Source Number-Counts 



By integrating the predicted flux, arriving from different directions, one can construct 
maps of the simulated 7-ray sky. Maps of the entire sky (Figure 8) and of two selected 
regions (Figure 9) are presented in the following. They reveal a globally isotropic picture, 
with strong local fluctuations, e.g. flux variation over 2 orders of magnitude at ~ 1° resolution 
and over more than 3 orders of magnitude at ~ 5' resolution. Such maps of the predicted 
7-ray sky can be directly compared to the maps produced by detectors such as EGRET 
(Sreekumar et al. 1998) aboard the CGRO satellite, detectors on board GLAST, planned 
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for launch in 2006, and Cherenkov telescopes such as MAGIC, VERITAS 6 and HESS 7 . 
Since such devices have finite angular resolution, the simulated maps must be convolved 
with appropriate filter functions, as demonstrated in Figure 9. The relevant parameters of 
EGRET, GLAST and MAGIC are summarized in Table 4. 

Various 7-ray sources, of different shapes and sizes, are observed in the simulated maps 
of the 7-ray sky. Some of these sources exhibit irregular shapes, revealing the underlying 
structure of their emitter. Two typical examples, depicted in Figure 9, are strong emission 
from a semi-spherical source, and weak emission from a filamentary object. Some interest- 
ing features of the semi-spherical 7-ray source are shown in Figure 10. This source is the 
signature of a nearby rich galaxy cluster, lying ~ 53 Mpc from the simulated observer, at 
a redshift of z ~ 0.012. The 7-ray emission from the cluster appears approximately as an 
elliptic ring, corresponding to a diameter 5 — 10 Mpc, surrounding the location of the cluster 
and indicating the position of its accretion shock. The luminosity has two evident peaks 
along the circumference of the ring, approximately at the locations of its intersections with 
two large galaxy filaments. This suggests that the strong emission from these two regions 
is due to shocked gas, channeled by the filaments into the cluster region. Other galaxy fila- 
ments may be responsible for less luminous peaks of emission along the circumference of the 
ring, on the right hand side of Figure 10 images. 

The simulated sky maps produced may be decomposed into discrete 7-ray sources, 
yielding number counts of sources above given photon energies (Figure 11). This is performed 
using a greedy source-identification algorithm devised for the purpose, discussed in Appendix 
B. The resulting number counts of sources above 100 MeV and above 1 GeV fall short of the 
EGRET detection threshold. Our simulation thus predicts that none of the ~ 60 unidentified 
extragalactic EGRET sources (Ozel & Thompson 1996) can be attributed to intergalactic 
shock-induced emission. The number-counts are especially sensitive to the fraction of shock- 
energy transferred to the relativistic electrons. For £ e = 0.05 our simulation predicts a few 
dozen sources detectable by GLAST at both photon energy ranges, among them 16 well- 
resolved sources above 1 GeV associated with different objects in the sky (see Figure 8). For 
£ e = 0.02, the simulation predicts only one source detectable by GLAST, whereas £ e = 0.10 
yields more than a hundred GLAST sources, well resolved above 1 GeV, yet too faint for 
EGRET detection. 

Since our source-identification algorithm was devised for spherical sources, very large, 
nearby objects tend to be broken down and identified as several separate, weaker ones. 



6 See http://veritas.sao.arisona.edu/ 

7 See http://mpi-hd.mpg.de/hfm/HESS/HESS.html 
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This implies that when comparing our source number counts with analytical predictions 
(Waxman & Loeb 2000; Totani & Kitayama 2000), one should impose a threshold on the 
angular size of the sources. As evident from Figure 11, Our number counts fall short of those 
of Waxman & Loeb (2000) by a factor of ~ 6. This may be explained, at least partially, 
by the same factors listed above that reduce the simulated diffuse 7-ray flux compared to 
the analytic calculation: the lower average temperature of the simulation suggests less hot, 
massive clusters, and combined with the thermal efficiency of the shocks and the fraction 
of electrons that had sufficient time to cool, both lower than assumed in the analytical 
predictions, may account for the discrepancy factor. 



6. Discussion 



We have studied the 7-ray emission from intergalactic shock waves, using a hydrody- 
namic ACDM cosmological simulation, according to the model proposed by Loeb & Waxman 
(2000). A simple approach for extracting large-scale shocks from a Lagrangian simulation, 
based on following entropy changes of single simulated particles, was employed. This method 
was shown to identify most of the strong (Mach number M > 100) accretion shocks in the 
simulation containing a large fraction (/ sh ~ 2/3) of the gas thermal-energy gain, although 
some of the weaker (M < 10) flow or merger shocks might be missed. Relativistic electron 
populations, Fermi-accelerated by the shocks, were injected into the shocked gas, assuming 
that a fraction £ e = 0.05 of the shock thermal energy was transferred into relativistic elec- 
trons. The inverse- Compton radiation, emitted as these electrons scatter CMB photons, was 
then calculated and integrated over space and time, to yield the radiation detected by an 
imaginary observer. 

This radiation was found to extend from the optical band and into the far 7-ray regime, 
containing roughly equal energy flux per decade in photon energy, e 2 [dJ 1 (e)/de 1 ] ~ 50 — 
160 eV cm -2 s" 1 sr -1 , in the photon energy range 10 keV < e < 1 TeV. The calculated 
spectral slope, dJ/de ~ t~ 2m for 10 keV < e < 1 GeV and dJ/de ~ e" 213 for 10 GeV < e < 
5 TeV (before accounting for photon-photon pair-production cascade), is consistent with the 
predictions of Loeb & Waxman (2000) and with the EGRET observations. The energy flux 
constitutes ~ 10% of the EGRB flux, or up to ~ 15% of the EGRB flux after subtracting 
the expected contribution from unresolved point sources based on empirical modeling of the 
luminosity function of blazars (Mukherjee & Chiang 1999). 

The simulated 7-ray sources, mostly spheroidal halos associated with large-scale struc- 
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ture, fall short of the EGRET detection threshold, thus producing a diffuse background. 
Hence, the simulation predicts that none of the ~ 60 unidentified EGRET extragalactic 7- 
ray sources can be attributed to intergalactic shock-induced radiation. Several sources do fall 
within the detection range of GLAST, planned for launch in 2006, provided that £ e > 0.03. 
The number of detectable sources is sensitive to the fraction of shock-energy transferred to 
the relativistic electrons: For £ e = 0.05 we find several dozen sources detectable by GLAST, 
among them 16 well-resolved sources (see Figures 8, 11), whereas £ e = 0.10 is still insufficient 
for EGRET detection, but leads to more than a hundred well-resolved GLAST 7-ray sources. 
Such sources could be identified as emission from intergalactic shock-accelerated electrons, 
by association with large-scale structure, a characteristic dJ/de ~ e -2 spectrum, extending 
up to the far 7-ray regime, and a low cutoff at photon energy e min ~ 5 (t/10 9 yr) keV, 
caused by low-energy electrons not having sufficient time to cool during the lifetime t of 
the source. Detection of intergalactic shock-induced 7-ray sources may be used to calibrate 
£ e , once their temperature has been determined. The shock-induced 7-ray emission, a first 
direct tracer of the warm-hot IGM (10 5 K < T < 10 7 K), may thus be used to determine 
the baryon density in the present-day universe. 

We have examined the 7-ray signature of a nearby simulated rich galaxy cluster, lying 
at a redshift z ~ 0.012. An elliptic emission ring of a diameter corresponding to 5 — 10 Mpc 
surrounds the cluster, tracing the location of its accretion shock (see Figures 9, 10). The 
luminosity peaks along the circumference of this ring, probably at the locations of intersec- 
tions with galaxy filaments, channeling gas into the cluster region. Similar features should 
be detected by GLAST or by the MAGIC telescope in nearby rich galaxy clusters, at pho- 
ton energies above 10 GeV, and by VERITAS and HESS above 200 GeV, if the background 
signal is assumed flat. A ring-like feature resembling the simulated emission should thus be 
detected in such galaxy clusters lying at redshifts z < 0.01, whereas the brightest regions of 
emission along the circumference of the ring could be detected in clusters with redshifts up 
to z ~ 0.025. These bright regions of emission lie above the background signal and should 
be detected at redshifts z < 0.015 even for £ e = 0.02, whereas £ e = 0.10 will enable the 
detection of the emission ring itself at similar redshifts. Detection of such a ring will enable 
a direct study of the cluster accretion shock, as well as establish the validity of our model. 
Detected bright regions along the circumference of such a ring may be tested to comply with 
data regarding a nearby galaxy filament and be used to measure the density and velocity of 
the untraced gas in such a filament. 

As mentioned above, our algorithm for extracting shocks succeeded in capturing strong 
shocks, but missed some of the weaker (M < 10) ones. Accounting for these weak shocks 
could only have a minute effect on the predicted diffuse 7-ray radiation, since most (~ 2/3) 
of the thermal energy has been identified in strong shocks and because the steep spectrum of 
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electrons accelerated by weak shocks reduces their 7-ray efficiency. Weak shocks could, how- 
ever, increase the 7-ray brightness of nearby structures such as galaxy clusters encompassing 
ongoing merger processes, mainly in low photon energies, slightly altering the appearance of 
their inner regions and enhancing the predicted source number counts. 

In our work we have neglected the radiation resulting from the protons accelerated in the 
intergalactic shocks. These protons lead to 7-ray emission as they scatter off thermal protons 
to produce pions, most of the 7-ray flux originating from 7r° decay. Even with a large proton- 
to-electron energy ratio 77, the radiation resulting from the relativistic protons is, on average, 
negligible compared to the emission from electrons. Only in the densest (n > 10~ 3 cm -3 ) re- 
gions in the cores of clusters, the relativistic protons may alter the 7-ray spectrum. This effect 
may make the cores of clusters even brighter 7-ray sources, provided that is not eliminated 
by diffusion of the relativistic protons out of the over-dense regions. For the values of £ e and 
77 discussed above, the relativistic protons are dynamically significant and play an important 
role in shock dynamics and in particle acceleration. These effects, possibly detected in the 
radio spectra of some supernovae remnants, are small (Reynolds & Ellison 1992), suggesting 
that the linear approach adopted in this paper is sufficiently adequate. 

One can confirm that at least part of the EGRB is associated with the large-scale struc- 
ture of the universe, by cross-correlating 7-ray maps with other sources that trace the same 
structure, such as galaxies, X-ray gas and the Sunyaev-Zel'dovich effect. The intergalactic 
shock-origin of 7-ray radiation may be verified by cross-correlating the EGRB with radio 
maps, partly attributed to synchrotron emission from the same shock-accelerated electrons. 
The synchrotron radio radiation is sub-dominant to the CMB, although its strong fluctua- 
tions on sub-degree scales were shown to contaminate CMB fluctuations at low [y < 10 GHz) 
frequencies (Waxman & Loeb 2000). Detection of a radio - 7-ray cross-correlation signal 
from nearby clusters may be used to measure the intergalactic magnetic field. The auto- 
correlation of inverse-Compton or synchrotron radiation, as well as their cross-correlation, 
may be calculated directly from our simulation. 

A study of the shock-induced component of the extragalactic radiation, combined with 
a calculation of shock-induced emission according to cosmological simulations as described in 
this paper, may provide an important test of structure-formation models. Angular statistics 
of the 7-ray radiation, both on its own and cross-correlated with other sources, may be 
measured and compared to the predictions of such simulations. This will enable a study 
of the large-scale structure of the universe, and provide insight to the structure-formation 
process itself. We plan to apply the tools, developed in this research, to other cosmological 
simulations, including high-resolution simulations incorporating radiative cooling processes 
and supernovae feedback (e.g. Springel & Hernquist 2002b). Such research will permit 
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an extensive study of shock-induced radiation, in the framework of various cosmological 
structure-formation models. Other physical phenomena, such as the acceleration of cosmic- 
rays by intergalactic shocks, may thus be investigated as well. 

Order-of-magnitude estimates, presented in §2.3.1, predict an inverse-Compton flux 
higher by a factor of ~ 10 than obtained from the simulation we have studied. This is 
primarily attributable to the lower present-day temperature reached by the simulation, 
T ~ 4 x 10 6 K, and to its slower heating-rate, which seem to comply with other recent 
simulations. Our simulation neglected radiative cooling, suffered from relatively low mass 
resolution (~ 1O 11 M ), and did not include feedback from supernovae explosions, photo- 
heating the surrounding matter. The effect of radiative cooling on the accretion shocks is 
expected to be small, as cooling is efficient in hot, dense regions, beyond these shocks, and 
in cold, dilute regions, where the temperatures are low to begin with. Cooling should have 
a non-negligible effect on merger and flow shocks, but since these are much weaker than 
the accretion shocks, this should not significantly alter our results. In order to examine the 
sensitivity of our results to the low resolution of the simulation, a preliminary check of a 
similar adiabatic ACDM simulation, with identical cosmological parameters but with eight 
times the mass resolution, was carried out, exhibiting no substantial difference from the 
results presented above. Outflows from galaxies (due to starburst or quasar activities) could 
also result in shock acceleration of relativistic electrons as they impact on the surrounding 
IGM. At the same time, preheating of the collapsing matter by, for example, supernovae, 
may lead to weaker, cooler shocks, as recently pointed out by Totani & Inoue (2001), and 
thus result in softer, weaker photon emission. However, Dave et al. (2001b) find that the 
luminosity-temperature relation for galaxy clusters, which originally motivated the addition 
of preheating to structure formation models, may in fact be attributed to radiative cooling. 
Therefore, preheating may not have played an important role in structure formation after 
all, and thus have little effect on our results. 

Note added in proof: Recently, Scharf & Mukherjee (2002) have pointed out a possible 
association of 7-ray emission with the locations of galaxy clusters, by cross-correlating high 
Galactic latitude EGRET data with Abell clusters. The correlation signal, identified at a 
confidence level of 3a, is broadly consistent with the extended emission and with the source 
number counts discussed above, suggesting good prospects for the detection of extended 
7-ray emission from clusters by future observations. 
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A. Integration Parameters 



In order to perform the space-time integration schemes described in section 4.3.2, one 
must find the redshift dependence of several parameters: the time that elapsed between 
redshift z and the present, the coordinate distance, traversed by a photon in that time, and 
the minimal Lorentz factor of electrons that lost a significant fraction of their energy since 
redshift z. Formulae for these parameters are provided below. 

Using the FRW equations, we find that the time that elapsed between redshift z and 
the present is given by: 



At = t - t(z) 



1 



R(z=0) j 

R(z) R dR ~ H ° J Q (1 + Z)h(z) 



dz , 



(Al) 



where h(z) = H(z)/H . The coordinate distance traversed by a photon in that period is 
given by: 
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where g(x) is a function of the curvature k, defined by: 
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For our ACDM model, where £7a + = 1 and h(z) = a/^a + (1 + z) 3 Qm, the above 
equations are solved to give: 
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where (5 is the incomplete beta function, defined by: 
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In a similar fashion, we find the redshift dependence of 7 min (z,p), the minimal Lorentz 
factor of an electron, that looses a fraction p of its energy between redshift z and today 
by inverse-Compton scattering. This parameter, which determines the low cutoff of the 
resulting spectrum, is given by: 



7min(^,p) 



P 



3m e cH 



(1 -p) Aa T u CMB (z = 0) 
which for our ACDM model can be solved to give: 
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where F(x) is the hyper-geometric function 2 F\ (|, |, |, — x), defined by 2-fi (§, |, |, — x) = 

The redshift dependence of the parameters calculated above, also demonstrating the 
criteria for snapshot selection (simulation-box light-crossing time), is depicted in Figure 12. 



B. Source Identification Algorithm 



In order to count discrete 7-ray sources, predicted by the simulation, we have devised 
a greedy algorithm for extracting point sources from a sky map. This algorithm is used to 
isolate such sources, determine their location and flux, and measure their angular size. The 
algorithm operates iteratively, identifying sources and removing them from the map, until 
the contrast of the remaining map is smaller than some imposed threshold (of order 2). At 
each iteration, the algorithm searches for the brightest point in the sky, and attempts to 
identify the brightest source associated with it. Our candidate for the source is assumed to 
have its kernel positioned about this bright point. We determine the area around the kernel, 
associated with the source by examining concentric rings about the kernel and adding them 
to the 'growing' source. This accumulation process lasts as long as the flux within such 
a ring is sufficiently higher than the background, throughout its circumference, and yet 
monotonically decreasing, at least at one point along the circumference. Once the entire 
source has been identified, the flux associated with it is carefully removed from the map. 
The flux subtracted from each ring is determined according to the faintest point along its 
circumference, in order to prevent upsetting adjacent sources. 
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The source number counts produced by the algorithm were analyzed in section 5.2, 
with respect to the experimental devices EGRET and GLAST. In order to do so, algorithm 
parameters were tuned to mimic the source identification process, used to analyze the ex- 
perimental data. Hence, the angular resolution of the analyzed sky maps was set to match 
the point-source resolution of the experimental devices, and the angular size of each source 
was compared to the photon-position error of each experiment, as illustrated in Figure 11. 

Figure 13 demonstrates the results of the algorithm, when applied to two regions in the 
simulated 7-ray sky. The sky maps displayed in this paper were produced by approximating 
the source, associated with each SPH particle, as a rectangle in — <fi space. The resulting 
image of a source, produced by a single SPH particle, thus assumes a trapezoidal-like shape 
when mapped to a surface, centered at the source. Similarly, The concentric 'rings' discussed 
above are in fact rectangles in 9 — space. The identified sources are thus depicted in Figure 
13 as ellipses instead of rectangles for illustrative purposes only. 
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Fig. 1. — Imaging a partial slice of the simulation box at z — 0. Width and height are 
100 Mpc, depth is 10 Mpc. 

left: Maximal baryon number density through the slice. Color scale: log 10 (n [cm" 3 ]). 
right: Maximal temperature through the slice. Color scale: log 10 (T [K]). 
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Fig. 2. — Thermodynamic evolution of the simulated universe at recent epochs. 
left: The temperature averaged over mass (p 1 averaging, solid line) and over volume (p° 
averaging, dashed line). The mass averaged temperature according to the approximation 
of equation (5) is plotted as well, for k — 1 (as assumed in §2.1, dash-dotted line) and for 
k = 0.67 (best fit to data, dotted line). 

right: Baryon number density averaged over mass (solid line), over volume (dashed line) and 
over temperature (dash-dotted line). The volume-averaged density overshoots its predicted 
value (~ [1 + z] 3 , dotted line), probably due to limited simulation resolution in dilute regions. 
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Fig. 3. — Radiative cooling of the simulated gas. 

left: Mass distribution at z = in the phase space of temperature and baryon number density. 
Gray-scale: the logarithm (base 10) of the mass fraction m~J:Am/ [Alog 10 (T) • Alog 10 (n)]. 
Superimposed black contours denote the cooling time t cool due to Bremsstrahlung and line 
emission (Peacock & Heavens 1990), normalized to the Hubble time t H . 
right: Mass fraction of particles with various normalized cooling times in different epochs. 
Gray scale: t cool /t H (z). 
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Fig. 4. — Statistics of normalized entropy-change, AS 1 *, accumulated by simulated particles. 
Plotted are the mass fraction f(AS'*) = m~J:dm(AS'*) / dAS* (thin solid line), integrated 
mass fraction m^ 1 m(AS t > AS'*) = f^L, f(AS*)dAS* (thick solid line) and the shock 
Mach number, M, corresponding to positive AS* values (thick dashed line, right axes). 
left: Total entropy accumulated in the recent epoch, AS* — S*(z = 0) — S*(z = 2). The dis- 
tribution of negative AS"* values, approximated by a normal distribution (thin dashed line) , 
is likely due to numerical errors, leaving behind a narrower distribution (well-deconvolved 
up to AS** ~ 1, dotted line). High positive AS* values are also well-fitted by a normal 
distribution (dash-dotted line). 

right: Average AS* between two consecutive snapshots. A minimal threshold for shock de- 
tection (dash-dotted line) was imposed at AS* = 0.4, where estimated mass fraction due to 
numerical noise is f(AS*) ~ 1% (according to negative AS* distribution). On average, 10% 
of the gas experiences higher entropy increases between consecutive snapshots. 
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Fig. 5. — Results of shock extraction between z = 2 and z = 0. 

left: Redshift dependence of the comoving energy densities (left axes), accumulated since 
z — 2, of: the simulated gas (thick solid line), the identified shocks (thin solid line), shock- 
accelerated electron population assuming £ e = 0.05 (thick dashed line), electrons within the 
relevant energy range j mia < 7 < 7 max (thin dashed line) and the redshifted inverse- Compton 
radiation reaching an observer today (dash-dotted line). Also shown (right axes) are the 
fraction of the mass that was shocked at least once since z = 2 (thin dotted line), and the 
fraction of the mass processed by shocks (n shocks experienced by the same particle counted 
as n times its mass, thick dotted line). 

right: Mass fraction f(AS*) with entropy change AS* accumulated between z = 2 and z — 0, 
for all SPH particles before shock extraction (dashed line), among the particles identified 
as shocked (solid line), and among these particles before multiple shocks experienced by 
the particles were separated (dotted line). Also shown (right axes) is the Mach number M 
corresponding to AS* (dash-dotted line). 
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Fig. 6. — The gas shocked between two consecutive snapshots, taken at z — 0.14 and at 
z = 0.09. Color scale shows the maximal baryon number density through a slice of the 
simulation box, in log^n [cm -3 ]). 

left: The same slice as shown in Figure 1, with dimensions 100 Mpc x 100 Mpc x 10 Mpc. 
right: Enlarged region in the left image, marked there by a bright square, with dimensions 
20 Mpc x 20 Mpc x 5 Mpc. 
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Fig. 7. — Inverse- Compton spectrum (solid line), contributing ~ 10% of the EGRET- 
measured flux (error bars, Sreekumar et al. 1998) and the estimates made in §2.3.1 (dashed 
line). Below 10 MeV the observed flux increases significantly: the shaded region represents 
results of the Solar Maximum Mission (SMM, Watanabe et al. 1997). The expected contribu- 
tion from unresolved point sources, based on empirical modeling of the luminosity function 
of blazars (Mukherjee & Chiang 1999), is shown (circle with error bar). Pair production 
cascade effectively cuts off the spectrum at photon energy e ~ 1 TeV (dotted line), slightly 
enhancing the spectrum below this energy. The approximated spectra, equation (57) and 
equation (58), are also shown (dash-dotted lines). The step-like feature at low energies is a 
consequence of the temporal resolution of selected snapshots (and corresponding j min values). 
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Fig. 8. — Full 7-ray sky map at photon energies above 100 MeV. Color scale represents 
log 10 [J/ J), where J = 1.1 x 10~ 6 cm -2 s _1 sr _1 is the average photon number flux at these 
energies. Longitudes and latitudes (dark lines) are plotted 45° apart. Bright contours enclose 
regions magnified in Figure 9. Blue circles show simulated sources that will be well-resolved 
by GLAST in photon energies above 1 GeV (for £ e = 0.05). Angular resolution is ~ 42'. 
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Fig. 9. — Photon number flux above 1 GeV from two 16° x 16° regions in the sky, marked 
in Figure 8 by dashed (left images) and solid (right images) bright contours. Images are 
convolved with Gaussian filter functions, simulating photon angular spread of standard de- 
viation 1.5° (for EGRET, upper images) and 0.42° (for GLAST, bottom images). Color scale 
represents log 10 (J/J), where J = 9.9 x 10~ 8 cm~ 2 s _1 sr _1 is the average photon number 
flux above 1 GeV. 
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Fig. 10. — Emission from a nearby rich galaxy cluster. 

upper, left: Photon number flux above 10 GeV from the 16° x 16° region shown in Figure 9 
(right images). The image is convolved with a Gaussian filter function of standard deviation 0.2°, 
the expected angular resolution of MAGIC in these energies. Color scale: log 10 (J/J), where 
J = 8.2 x 10~ 9 cm -2 s _1 sr _1 is the average photon number flux above 10 GeV. 
upper, right: Maximal baryon number density through a thick slice of the simulation box at z = 0, 
showing a rich galaxy cluster and connected galaxy filaments. The dashed frame marks the central 
~ 15 Mpc x 15 Mpc region, the emission from which is shown in the upper left image. The slice, of 
dimensions 32 Mpc x 32 Mpc x 8 Mpc, is oriented perpendicular to the line of sight, with its center 
located 53 Mpc {z ~ 0.012) from the observer. Color scale: log 10 (n [cm -3 ]). 

bottom: Maximal temperature through the slice. Color scale: log 10 (T[K]). Contours trace lines of 
constant number flux above 10 GeV, higher than the average flux in this energy range by a factor 
ranging from 1.7 (black lines) to 22 (green lines). 
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Fig. 11. — Cumulative all-sky number of 7-ray sources with observed photon number 
flux exceeding F, for photon energies above 100 MeV (left panel) and above 1 GeV (right 
panel). Plotted are the number of all sources (thin solid line) and the number of sources 
with angular size smaller than the la photon-arrival spread of EGRET (dash-dotted line) 
and of GLAST (thin dashed line). Estimates based on the Press-Schechter mass func- 
tion (Waxman & Loeb 2000) predict a higher number of bright sources (thick solid line) 
and bright sources with angular size < 1° (thick dashed line). In both energy ranges, 
the number of the brightest sources as a function of flux is well-approximated by a power 
law iV(> F) ~ F~ 2 (open circles), steeper than the Press-Schechter based estimates 
N(> F) ~ F- L38 . Also plotted are the detection flux thresholds of EGRET and GLAST 
(vertical dashed lines). 
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Fig. 12. — Redshift dependence of the parameters, required for EGRB integration. 
left: Coordinate distance (rR , thin solid line) between a source emitting a photon at redshift 
z and an observer detecting the photon today (in 10 10 ly units), the angular diameter dis- 
tance (ri? / [1 + z], dashed line) and the luminosity distance (rR [1 + z], dash-dotted line) 
between them. The time that elapsed between emission and detection is given in 10 10 yr 
units (thick solid line). Circles and dotted lines illustrate the choice of simulation snapshots 
examined, separated by the comoving simulation-box light-crossing time. 
right: Minimal Lorentz factor of electrons, emitting at least half their energy between red- 
shift z and the present (dashed line), and maximal Lorentz factor reached by the shocked gas 
according to the simulation (solid line) and according to equation (12) (dash-dotted line), 
where a post-shock magnetic field of magnitude B = 0.1 fiG was assumed. 






Fig. 13. — Results of the source-identification algorithm, when applied to two regions in the 
7-ray sky. Both images are of angular size 16° x 16°, and are similar (but not identical) to 
those presented in Figure 9. Color scale is identical to the one shown in Figure 9. 
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Table 1. Parameters of the cosmological model and linear structure formation theory. 



Parameter 


71 /T 

Meaning 


t r t 

Value 


h 


Hubble parameter 


0.67 


k 


Curvature 







Matter energy density 


0.3 




Dark matter energy density 


0.26 




Baryon energy density 


0.04 


n A 


Vacuum energy density 


0.7 


X 


Hydrogen mass fraction 


0.76 


n 


Fluctuation spectrum slope 


1 




Spectrum normalization 


0.9 


^start 


Starting redshift for the simulation 


50 



Table 2. Code parameters of the simulation. 



Parameter 


Meaning 


Value 


L coul R{t ) a 


Simulation box comoving length 


134/1" 1 Mpc ~ 200 Mpc 




Number of dark matter particles 


224 3 ~ 11 x 10 6 


M DM 


Mass of a dark matter particle 


2.3 x 10 10 M 


N SPH 


Number of gas (SPH) particles 


224 3 ~ 11 x 10 6 


M SPH 


Mass of an SPH particle 


3.55 x 10 9 M 


^min 


Minimal SPH softening 


Qh^ 1 kpc ~ 9 kpc 


N s 


Number of SPH neighbors in kernel 


32 


M ies 


Minimal Mass Resolution 


~ 10 11 M 


S G R(t ) a 


Gravitational softening 


25/1- 1 kpc ~ 37 kpc 



a i?(t ) is the present-day value of the cosmological scale factor. 
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Table 3. Comoving energy density u com of various simulated components. 



Parameter 


Component 


M co ™ [eV cm 




Ui 


gas, at z = 2 


2.77 x 10" 


-5 


u f 


gas, at z = 


1.72 x 10- 


-4 


Au = Uf — Ui 


gain of gas during < z < 2 


1.44 x 10- 


-4 




detected shocks 


9.62 x 10- 


-5 


u e 


accelerated electrons (£ e = 0.05) 


4.81 x 10- 


-6 


^eVTmin *^ 1e *^ Tmax) 


electrons with 7 min < 7 e < 7 max 


2.34 x 10- 


-6 


M 7 


radiation, redshifted 


1.40 x 10- 


-6 



Table 4. Parameters of EGRET and of the conceptual designs of GLAST and of MAGIC. 



Parameter 


Photon Energy 


EGRET a 


GLAST a 


MAGIC b 






(one year all-sky) 


(one year all-sky) 


(50 hours) 


Point source 


E > 100 MeV 


5.4 x 10~ 8 


1.5 x 10~ 9 




- sensitivity 


E > 1 GeV 


1.2 x 10~ 8 


1.5 x 10- 10 




[cm- 2 s _1 ] 


E > 10 GeV 




1.0 x 10~ 10 


1.0 x 10~ 10 


Photon - 


E > 100 MeV 


5.6° 


2.5° 




position - 


E > 1 GeV 


1.5° 


0.42° 




error (la) 


E > 10 GeV 




0.1° 


0.2° 



a See http://www-glast.stanford.edu 

b See Gonzalez et al. (1997) and http://hegral.mppmu.mpg.de/MAGICWeb 



